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Abstract 

We present a practical algorithm based on symplectic splitting meth¬ 
ods to integrate numerically in time the Schrodinger equation. When 
discretized in space, the Schrodinger equation can be recast as a classi¬ 
cal Hamiltonian system corresponding to a generalized high-dimensional 
separable harmonic oscillator. The particular structure of this system 
combined with previously obtained stability and error analyses allows us 
to construct a set of highly efficient symplectic integrators with sharp er¬ 
ror bounds and optimized for different tolerances and time integration 
intervals. They can be considered, in this setting, as polynomial approx¬ 
imations to the matrix exponential in a similar way as methods based on 
Chebyshev and Taylor polynomials. The theoretical analysis, supported 
by numerical experiments, indicates that the new methods are more effi¬ 
cient than schemes based on Chebyshev polynomials for all tolerances and 
time intervals. The algorithm we present incorporates the new splitting 
methods and automatically selects the most efficient scheme given a tol¬ 
erance, a time integration interval and an estimate on the spectral radius 
of the Hamiltonian. 
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1 Introduction 

When investigating the dynamical behavior of quantum systems of low to mod¬ 
erate dimension, very often it is necessary to solve numerically the time depen¬ 
dent Schrodinger equation (h = 1) 

= Hip(x,t), q(x,o) = VAA (1) 
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Here H is the Hamiltonian operator, if; : x M —> C is the wave function 

representing the state of the system and ipo(x) is the initial state. For simplicity, 
in the sequel we consider H = T + V, with the kinetic energy operator T = 
— A/(2/z) for a reduced mass fj, > 0 and a potential V. although the procedure 
presented in this paper is also valid for more general Hamiltonian operators. 

The solution of (jTj) can be expressed as 

fp(x,t) = U(t)ip 0 (x), (2) 

the (unitary) evolution operator U being formally given by U(t) = e~ ltH . In 
practice, however, it is not possible to get a closed expression for U(t), and 
so numerical methods are applied to get reliable approximations. This process 
involves typically two stages. In the first a discrete spatial representation of 
the initial wave function ipo(x) and the operator H on an appropriate grid are 
constructed. In the second, this finite representation is propagated in time with 
a numerical integrator. 

As for the space discretization process, several techniques can be used, 
depending on the particular problem one aims to analyze: finite difference 
schemes, spectral methods based on collocation with trigonometric polynomi¬ 
als, Galerkin method with a Hermite basis, etc, both in one or more dimensions 
(see m and references therein). The space discretization process restricts the 
energy range of the approximation and imposes an upper bound to the high 
frequency components represented by the discrete solution. 

In any event, once this process has been carried out, one has the linear 
system of ordinary differential equations 

i—u(t) = Hu(t ), u(0) = uq £ C N , (3) 

where u(t) now represents a discretized version of the wave function ip(x, t ) at 
the N space grid points, with N usually a large number. The goal is then to 
compute u(t) at a given target time t from the known value of u( 0) = uq. The 
N x N matrix H (and in particular its discrete spectrum) depends of course on 
the particular space discretization carried out. We will hereafter assume that 
H is a real symmetric matrix which implies that it can be diagonalized with 
real eigenvalues. 

The exact solution of eq. ([3]) reads 

u(t) = e~ ltH u 0 , ( 4 ) 

but computing the matrix exponential e~ ltH by diagonalizing H (usually, a ma¬ 
trix of large dimension and large norm) is prohibitively expensive. An effective 
alternative consists in computing approximations of u(t) of the form 

u(t) « P m (tH)uo, (5) 

where P m (y) is a polynomial in y that approximates the exponential e~ iy , since 
in that case only multiplications of the matrix H with vectors u are necessary. 
These products can be efficiently evaluated in complex variables (provided that 
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a Fourier spectral method is used to obtain the discretized version ([3]) of Q) 
with the complex-to-complex Fast Fourier Transform (FFT) algorithm [5J flT| 

nil ns]. 

There are different choices for such a polynomial P m (y). For instance, one 
may consider truncated Taylor or Chebyshev series expansion of e~ zy for an 
appropriate real interval of y, or the Lanczos method, where the polynomial 
is determined by a Galerkin approximation on the Krylov space spanned by 
u 0 , Hu 0 ,H m ~ 1 u o EH], 

In this paper we consider yet another kind of polynomial approximation 
to e~ ltH uo, namely one based on explicit symplectic splitting methods ciia 
Huai. This approach can be applied under the same assumptions than the 
Chebyshev method, the main difference being the following. Whereas in the 
Chebyshev (or Taylor) method the approximation © is constructed by evaluat¬ 
ing products of the form Hu, where u E C^, with symplectic splitting methods 
one writes u = q + ip, q,p E . The algorithm then proceeds by successively 
computing real matrix-vector products Hq and Hp with different weights, so 
that the real and imaginary parts of e~ ltH UQ are approximated in a different 
way, with a much reduced computational cost. 

More specifically, if a spatial discretization based on Fourier spectral meth¬ 
ods is considered, then the cost of computing Hu, u E C , amounts essentially 
to one complex-to-complex FFT and its inverse, whereas in the case of Hv, 
v E M. N , one has to evaluate one real-to-complex FFT and its inverse complex- 
to-real FFT, and this process requires half the computing time of the fully 
complex case. As a result, the proposed algorithm based on splitting methods 
turns out to be between 1.4 and 2 times faster than the Chebyshev method for 
the same accuracy in all the examples we have analyzed. Moreover, the proce¬ 
dure is easy to implement and the resulting approximations preserve important 
qualitative properties of the exact solution. 

The algorithm we present here has embedded several symplectic splitting 
schemes designed according to different optimization criteria with the purpose 
of covering most of the cases one finds in practical applications (high accuracy 
over long time intervals, low accuracy over short times, etc.). The computation 
of the coefficients of the methods, which constitutes a non-trivial task by itself, is 
largely based on the stability and error analysis of splitting methods carried out 
in m i- Given a target value of time t and an error tolerance, the algorithm 
selects a specific symplectic splitting scheme leading to a numerical solution 
with the prescribed accuracy and the minimum computational work, measured 
as the number of real matrix-vector products. By construction, the algorithm 
developed here is aimed to be applied for the same problems and under the same 
assumptions as the Chebyshev method, with a remarkable gain in efficiency for 
all the examples we have tested. 

The plan of the paper is the following. Since our procedure may be con¬ 
sidered as an alternative to the Chebyshev method, in section [2] we summarize 
the main features of the schemes based on this polynomial approximation of 
In section [ 3 ] we analyze the stability and the global error of symplec¬ 
tic splitting methods in this context, and the actual algorithm is presented, 
whereas the comparison with Chebyshev (and Taylor as a reference) is carried 
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out in section [4] on a pair of selected numerical examples. 


2 Polynomial approximations 


2.1 General considerations 

Given a mth degree polynomial P m (y) approximating e ~ iy , the solution u(t) = 
e~ ltH uo of ([ 3 ]) at a prescribed target time t can be approximated as 

u(t) « ui = P m (t H) uq, (6) 


with the corresponding error (in Euclidean norm) bounded as 

\\ui — e~ ltH uo\\ < max \P m (t EA — e~ itEj \ \\uq\\ 

in terms of the (real) eigenvalues Eq,...,En-i of H. Assuming that the 
spectrum cr(H) = {Eq, ..., En-i} is contained in an interval of the form 
[A'min j l^max] j then 

IK - e~ ltH u 0 \\ < sup \P m (y) - e~ iy \\\u 0 \\. 

tE min <y<t ^max 


There are several possibilities to estimate £ maj( and E m ; n for different classes 
of matrices (see e.g. HU 2J. [22]). If PI can be decomposed as the sum 
H = T + V of two symmetric matrices with known lower and upper bounds for 
their eigenvalues, E m j n (resp. F max ) can be simply obtained as the sum of the 
lower (resp. upper) bounds of the eigenvalues of T and V. This happens, in 
particular, when the Hamiltonian operator H = —A/(2y,) + V is discretized by 
spectral Fourier collocation with N Fourier modes, in which case 

1 N 2 

Em in = min V (x) , F max = -— + max V (x) . (7) 

X Zfi 4 X 

In any event, once E m in and F max have been determined, we introduce 


E n 


+ E n 


a = 




Emax E m [ n 


H = H -al, 


( 8 ) 


so that the spectrum of the shifted operator H is contained in an interval 
centered at the origin, a(H) = {Eq — a ,..., F/v_ 1 — a) C [—/3,/3\. We thus 
have 


e~ itH 


Uq 


e ~ita e -iiH 


Uq. 


(9) 


Hence, we will hereafter assume without loss of generality that our problem 
consists in approximating e~ ltH UQ for a real symmetric matrix H with cr(H) C 
[— f3,/3\. In that case, 


\u 1 - 


e ltH UQ\ 


< e m (/3t) ||wo| 


( 10 ) 


where 

e m {0) = sup \e~ ty - P m (y)\. (11) 

-e<y<e 
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2.2 Taylor polynomial approximation 

The mth degree Taylor polynomial P^iv) corresponding to e~ iy is of course 


m (_A\k 

Pm(y) = E ~^y k 


k=0 


k\ 


( 12 ) 


and Horner’s algorithm provides an efficient way to compute u± 
namely 

Vo = uo 
do k = 1, m 

t 

Vk = u 0 i ■ - -Hyk—i 
m + 1 — k 

enddo 

ui = y m . 


P m( tH )u 0 , 

(13) 


The process requires storing three complex vectors (or equivalently, 6 real vec¬ 
tors). 

An error estimate of the form (10) can be obtained with e m (8) in @> 
replaced by its upper bound 




Qm -\-1 


('m + 1 )! ’ 


(14) 


Since ml 


\f2irm (rae) m for large values of m m , we can write 




:v/2 


7 xm 


/ 9e\ m+1 

J 


In consequence, we cannot expect to have a reasonably accurate approximation 
P^(tH)u o of e~ itH uo unless 


m > ed = e(31. 


In other words, increasing the value of the target time t where the solution is 
to be found and/or refining the spatial discretization (so that /3 gets larger) 
requires evaluating a higher degree Taylor polynomial. 


2.3 Chebyshev polynomial approximation 

The Chebyshev polynomial expansion scheme, proposed for the first time in 
the context of the Schrodinger equation in |19| . constitutes a standard tool to 
compute Q. A detailed analysis of the procedure, including error estimates for 
the problem at hand, can be found in [15]. For completeness, we review here 
some of its main features. 

The mth degree truncation of the Chebyshev series expansion of e - * y in the 
interval y G [—0, 6} is given by 


P%,e(v) = M8) + 2 Y^(-i) k MO)T k (y/ 0 ), (15) 

k=\ 
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Figure 1: Comparison of the required minimum polynomial degree m as func¬ 
tion of 0 = fit for Taylor (dashed line) and Chebyshev (continuous line) for 
different values of error tolerance: tol = 1(T 4 ,2 x 10 7 ,10 n . Diamonds, 
squares and circles stand for the computational cost (equivalent to a polyno¬ 
mial approximation of degree m ) for error tolerances below ICC 4 , 2 x ICO' and 
ICO 11 , respectively, obtained with symplectic splitting schemes in Table [lj 


where for each k, J k {t) is the Bessel function of the first kind |17j and T k {x) is 
the kth Chebyshev polynomial generated from the recursion 


T k+1 (x) = 2 xT k (x) -T k _ i(x), k > 1 (16) 


and To(x) = 1, T\(x) = x. According with the analysis in [15], e~ ltH uo can be 
approximated by p t (t H)uq with an error estimate of the form (10), where 
e m (0) in (JTTj) is replaced by its upper bound 


(Q\ _ 


( 0 ) = 4 


s l-0 2 /(2m+2) 2 


2m + 2 


m+1 


(17) 


In Figure [T] we depict the minimum degree m as a function of 6 = fit of 
Chebyshev approximations for prescribed tolerances tol = 10~ 4 , 2xl0~ 7 ,10 —11 , 
so that e())(/3t) < tol (continuous lines) in comparison with the corresponding 
degree m for Taylor approximations (dashed lines) such that e^[(/3t) < tol. 
Notice that Chebyshev always gives a similar accuracy with a lower degree 
polynomial (hence, with less computational cost), with a gain in efficiency of 
up to a factor of two for sufficiently large values of 6 = fit. 

Once the degree of the polynomial m has been chosen, given a certain error 
tolerance, target time t, and bound f3 of cr(H), one has to compute P mM tH ^ U 0 
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in an as efficient as possible way. This can be done with the Clenshaw recursive 
algorithm as follows: first evaluate the coefficients Ck = (—1 ) k Jk((3t) for k = 
0,1 ,,m and then compute recursively 

dm+2 — 0 , d m+ i — 0 

do k = m, m — 1,..., 1,0 

dk — Cfc no T ~pHdk -\-1 dk +2 (18) 

enddo 

ui = d 0 - d 2 , 

which produces u\ = p t (tH)u,o ~ e~ ltH uo as output. Clenshaw algorithm 
keeps only four complex vectors in memorjQ but the whole procedure has to 
be carried out for each value of m. Since the coefficients Ck are relatively small 
as k grows, the Clenshaw algorithm is stable and so it is possible to work with 
polynomials of very high degree (even in the thousands) provided the Bessel 
functions are accurately computed. 


3 Symplectic splitting methods 


3.1 General considerations 

An alternative to Chebyshev polynomial approximations of e~ ltH uo first con¬ 
sidered in 01 consists in applying specially designed splitting methods to 
numerically integrate the system Q recast in a more suitable form. 

By considering q = Re(u) £ M> N and p = Im(rt) £ M. N , equation ^ is 
equivalent to 

= (A + B)z, z(0) = z 0 , (19) 

where 



The solution z{t) = e t ^ A+B ' ) zo of (19) can be written in terms of the orthogonal 
and symplectic matrix 



cos (y) 
~ sin (y) 


sin (y) \ 
cos (y) ) 


( 21 ) 


as z(t) = 0(tH)zQ. To introduce general symplectic splitting methods in this 
setting, let us first show how the well known Strang splitting can be used to 
approximate e~ ltH UQ. Let m be a sufficiently large positive integer, so that for 
r = t/m , we consider the approximation 

e r(A+B) _ 

1 If the vectors are written in their real and imaginary part, and the algorithm is carried 
out in real variables, then the algorithm needs to store only seven real vectors instead of eight. 
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It is then clear that 


e^+B) = (V^) m « ( e i A e^ei A ) m = ei A (e" B 

or equivalently, 

0(tH) = e^ A+B) « K(tH) = e t°m+iA e tbmB e ta m A .. 
with 


e^) m 'i e rB e^ A , 

• e tftl B e tai A , (22) 


(fll, b \, Cl2) ■ • ■ > 0"m.i bm , ®m+l) — 


( 1 

1 

i 

i i i 

V 2m’ 

m ’ 

m ’ 

’ m ’ m ’ 2m 


(23) 


Due to the nilpotent structure of the matrices A and B in (20), the exponentials 

(24) 


in the definition (22) of K(tH ) take a particularly simple form, namely 

/ 1 


e * a j A = 


CLj t H 


I 


e t bj B _ 


-bjtH 


This analysis shows that the approximation K(tH)zo ~ e t ^ A+B ^Zo can be com¬ 
puted with the following procedure, similar in nature and equivalent in com¬ 
puting time to the Horner (13) and Clenshaw (18) algorithms: Given uq G C^, 


q := Re(u 0 ), 

p := Im(iio), 

do k = 1, m 

q := q + a k tHp 

p := p-b k tH q 

enddo 

q := q + a m+ it H p 
ui := q + ip, 


(25) 


producing v,\ ~ e~ ltH uo as output. Notice that it only requires storing three 
real vectors of dimension N (namely q, p, and w = Hp or w = Hq) instead 
of seven real vectors for the Clenshaw algorithm and six real vectors for the 
Horner algorithm. It is worth remarking that, since e A and e B are symplectic 
matrices, K{tH) is also symplectic. Unitarity is no longer preserved by this 
scheme, but neither the average error in energy nor the norm of the solution 
increases with time, since it is conjugate to a unitary method [2j. 

In practice, and in the same way as other polynomial approximations, it is 
convenient to apply Algorithm (25) with the original H replaced by the shifted 
version H considered in 


(and then make use of the equality (|9j)), so that the 
spectrum of H is contained in an interval of the form [— j3, f$] with f3 as sharp 
as possible. Therefore, in what follows we always assume that c(H) c [—/?,/?]. 

Although Algoritm (25) with coefficients (23) can be used in principle to 
approximate e~ ltH UQ, we next show that, for given values of m and 6 = /3t, 
much better approximations can be obtained if other sequences of coefficients 
(ai, bi, ci 2 , ■ ■ ■, a m , b m , a m+ 1 ) are chosen instead. To see how this can be done, 
an error estimate of the corresponding approximation (22) is necessary first. 











3.2 Error analysis 

For a given finite sequence of real numbers 


(«i, bi, «2 ’)••••> tirm ®m+1 ) ; 


( 26 ) 


Algorithm (25) produces an approximation of the form 


qi )=K(tH) qo 
Pi \Po 


(or equivalently, q\ + ip\ « e ltH (qo + ipo )) with 


(A+B) ( 9o 
Po 


K(tH) = 


K n (tH) K 12 (tH) 
K 21 (tH) K 22 (tH) 


(27) 


Here K\\ (y). K 22 (y ) are even polynomials of degree 2m, Ki 2 (y) and K 2 \{y) 
are odd polynomials of degree 2m — 1 and 2m + 1 respectively, and det K(y) = 
K u (y)K 22 (y) — K\ 2 (y)K 2 \(y) = 1. It is important to remark that for a given 


positive integer m, compared to Horner’s (13) and Clenshaw’s (18) algorithms, 


the degree of the polynomials involved in an m-stage splitting method (26) is 


twice the degree of the corresponding Taylor and Chebyshev polynomials, with 
the same computational cost. 

3.2.1 Error estimates for a single application of a splitting method 

We next focus on obtaining upper bounds for the error 

||(9i +ipi)-e- itH (q 0 + ip 0 )\\ = K(tH ) (^-0(tH) 

< \\K{tH)-0{tH)\\\\qo + ipo\\ 

in Euclidean norm. Since H is assumed to be a real symmetric matrix, it can 
be diagonalized as 


H = P 


( Eq 0 

0 E\ 


0 \ 
0 


0 0 0 

V 0 • • • 0 En-iJ 


p. 


where P is an orthogonal N x N matrix. We thus have 

K(t H) — 0(t H) = P t £P, 

where £ is the block-diagonal matrix (with 2x2 matrices at the diagonal) 


/ K(t Eq) — 0(t Eq) 0 

0 K{tEi)-0(tEi) 


0 

0 


\ 


0 


0 K(tE N _i) - 0(t.E N -i)J 
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and therefore 

\\K(tH)-0(tH)\\<\\S\\ = . max \\K(tEj) - 0(tEj)\\. 

Since \Ej\ < /3, j = 0,1,..., N — 1, we finally arrive at 

\\(qi + ipi) ~ e~ ltH {q 0 + ipo)\\ < e(fit) \\qo + ipo\\, (28) 

where 

e(9)= sup \\K(y) — 0(y)\\. (29) 

-9<y<6 

By taking into account that det K(y) = 1, the 2-norm of the 2x2 matrix 
K(y ) — 0(y ) can be explicitly computed to give 

\\I<(y) - 0(y)\\ = \/ ( c (v) ~ cos(y)) 2 + {S(y) - sin(y)) 2 

+ v / C(y) 2 + 5(y) 2 - 1, 

where 

C(y) = ^{Kn(y) + K 22 (y)), S(y) = ^( K 12 {y) - K 2l {y)). (30) 

Notice that det K{y) = 1 implies 

C{y) 2 + S(y) 2 - 1 = \(K u (y) - K 22 {y)f + \(K 12 {y) + K 2l {y)) 2 
and thus C(y) 2 + S(y) 2 — 1 > 0 for all real values of y. 


3.2.2 Error estimates for several steps of a splitting method 


Ideally, given a positive integer m and 9 = fit > 0, one would like to determine 
a sequence (26) of real numbers so that e(9) is minimized. The error bound 
e(9) being small implies that the (2m)th degree polynomial C(y ) (resp. the 
(2 m + l)th degree polynomial S(y)) is a good polynomial approximation of 
cos (y) (resp. sin (y)) for y £ [—9, 0], which implies that increasingly large values 
of 9 = (3t will require longer sequences of coefficients (that is, larger values 
of m), and consequently more computational work. The situation here is in 
complete analogy with what happened to Taylor and Chebyshev polynomial 
approximations in the previous section. 

By applying the methodology exposed in [3] we have determined several 


sequences (26) of length 2m + 1 of (near-to-optimal) coefficients for m up to 


60. The procedure is described in detail in the Appendix. As shown there, the 
task is by no means trivial, and severe technical difficulties arise when trying to 
extend the procedure to arbitrarily large values of 9 = fit (and hence arbitrarily 
long sequences of coefficients). This is in contrast with Taylor and Chebyshev 
approximations. 

This drawback can always be circumvented by approximating the solution 


z{t) = 0(t H)zq of the system of ordinary differential equations (19) in the 
standard step-by-step way. In our case, approximating z(t) in n steps of length 


r = 


n 


10 







simply consists in approximating 0(tH)zo = 0{nr H)zq = 0(rH) n ZQ by the 
vector K(t H) n zo, where K (y) is a 2 x 2 matrix with polynomial entries (defined 


in terms of the sequence (26) as before) that should approximate the rotation 
matrix 0(y) for y G [- yy, yy] ■ 

Clearly, the resulting procedure for approximating e~ itH uo can be written 


as an algorithm of the form (25), corresponding to a sequence of coefficients 


(with a (2m)-periodic pattern) of length 2 nm + 1. The corresponding error can 
be estimated as 

||(gi +ipi) -e~ ltH (q 0 + ip 0 )\\ < \\K(rH) n - 0{nrH)\\ ||g 0 + ipo\\ 

< e (n) (T/3)\\q 0 + ip 0 \\, (31) 


e v "'(9) = sup \\K(y) n — 0(ni 


where 

-8<y<e 

Our goal is then to minimize e^ n \6). A reasonable requirement is that K(y) n 
be bounded for all n. This only happens in general for a certain range of values 
of y. One thus defines the stability threshold y* as the largest non negative real 
number such that K(y) n is bounded independently of n > 1 for all y G (—y*, y*) 
[2], In particular, for the sequence (23) corresponding to the application of m 


steps of the Strang splitting, the stability threshold is y* = 2m. As a matter of 
fact, 2m is precisely the maximal stability threshold a sequence of coefficients 


(26) of length 2m + 1 can achieve m- 

From the analysis carried out in [3j, it is possible to show that 


\\K(y) n — 0(n\ 


< 2 sin(n(arccos(C'(y)) — y)/2) 

S(y? 


S{y) 2 l 

V 1 -C(y) 2 2 


-1 


.l-C(y ) 2 

provided that y G [—y*, y*]. This implies that, if t/3 < y*, then 

\\K{rH) n - 0{nrH)\\ < sup \\K(y) n - 0(ny)\\ = e (n )(r/3) 

—r/3<y<r/3 

< nn{TP) + (32) 


where 


y{9) = 

sup 


-8<y<G 

u(9) = 

sup 


- 8<y<e 


S(y)‘ 


l 

~ 1 + 2 


S(y) 2 


-1 


(33) 

(34) 


l-C(y) 2 

As mentioned before, we have determined several optimized splitting meth¬ 


ods of m stages (determined by a sequence of coefficients (26) of length 2m + 1 ) 


for m up to 60. The relevant parameters of such splitting methods are collected 
in Table 111 In this table, refers to a method of m stages, with error co- 

(1 3 ) 

efficients e(9), y(9), v(9) optimized for 9 = 7 m. For instance, method Mg 0 ' 
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can be used to approximate e~ ltH uo with an error bounded (according to (28) 
and Table[l]) by 1.2 x 10“ 9 ||rto|| provided that \t\ < 78/(3. Furthermore, e“ rf Mxo 
can be approximated by applying n steps of length t = t/n < r max := 78/(3 of 

f I of 

method Mg 0 ' ; with an error bounded (according to (32)) by 


(7.8 n x 10“ n + 1.2 x 10“ 9 


In some cases two methods with the same values of m and 7 = 0/m have 
been collected, in which case they are labeled a and b. For instance, methods 
Mg 4)a and Mgg’ 4 ^ 6 are both designed to approximate e~ ltH uo with n steps of 
length t = t/n < r max := 84/(3 of the method. However, they differ in the 


actual error estimate (32): in the first case, the error is bounded (provided that 


(3\t\ < 84 n) by (2.4n x 10 8 + 7.4 x 10 8 )||uo||, while the second one admits the 
error estimate (3.7 nx 10“ 9 + 2.6x 10 _6 )||ito||. This means that method A-lg 4 ^' 
will be more efficient if (3\t\ < 10452, and the opposite otherwise. 

Thus, given the upper bound (3 of the spectral radius of H and the target 
time t, if one wants to approximate e~ ltH uo by applying n steps of method 
Mg A)a 1 one should choose the smallest positive integer n such that 


t _ 84 

_ ^max •— > 

n p 


that is, n = Ceiling[f/3/84]. 


For instance, suppose the target time t and the bound (3 are such that t(3 = 1000. 
Then, clearly, n = 12, so that 12 steps of scheme Mg^ a have to be applied 
with step size r = 1000/(12 (3) — 83.3 /(3 to achieve the target time. In this way 
one gets an approximation with estimated error of size 3.62 x 10“' \\uo\\ with a 
computational work (2 nm = 2 x 12 x 60 = 1440 real matrix-vector products of 
the form Hv) comparable to the use of a Chebyshev polynomial approximation 
of degree 720. In contrast, to guarantee a similar precision with Chebyshev, a 
polynomial of degree at least 1135 is required, since this is the minimum value 
of m such that e,^(1000) ||uo|| < 3.62 x 10“' ||r£o||> with (-//(d) given in (JTt]) . 

It is worth remarking the error coefficients for Strang splitting method ( |23| ) 
with the same value of 7 = 9/m = 1.4 (also collected in Table [I]) are much 
larger than for methods Mgg' 4 ^ and Alg 4)h . 


3.2.3 Error estimates for combined splitting methods 

Sometimes it is just more efficient to apply a combination of two different 
methods instead of n steps of the same scheme. For instance, suppose that 
t(3 = 177 and we have an error tolerance of tol=10“'. If we use M^ { p a 
then t(3/ 84 ~ 2.1 so that, according with the previous considerations, method 
Mg q' 4) “ has to be used with n = 3 steps of size r = 59 /(3, much smaller than 
the value r max = 84/(3 for which the scheme has been designed. This would 
result in an approximation fulfilling the required error tolerance obtained with 
360 real matrix-vector products of the form Hv. A better strategy would be 
the following: apply two steps of scheme Mg q' 4 ^“ with step size r max = 84/(3 to 
approximate w = e“* 2Tmax h uq and then approximating e“* h~ 2rmax ) H w by using 
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some other method with less stages. More generally, we take n = Floor[f/3/84] 
steps of length r max = 84//J to get w = e~ inT[na * H uo and then we approximate 
e ~i (t-nTmax) h w w j^] 1 another method of Table [l] involving less stages. 

To decide which method has to be used for this last step, we need an er¬ 
ror estimate for the approximation obtained with such a combination of two 
methods. Assume that we apply n steps of length f of a method characterized 
by a 2 x 2 matrix K(y ) with polynomial entries, followed by a step of length 
r of a method characterized by the matrix K(y), where t = nr + r. From the 
preceding considerations, it is enough to estimate || K(t H)K(t H) n — 0(t H) ||. 
This can be done in terms of the functions fi(9), 9(9) associated to K(y) as 


defined in subsection 3.2.2, and the error function e(9) associated to K(9) as in 


subsection 3.2.1, together with the following function associated to K(y ): 


5{6) = sup \\K(y)\\-l. 

-0<y<6 


(35) 


Indeed, one obtains the following error estimate: 

\\K(t H)K(t H) n - 0(t H)\\ < \\K(t H) — e~ irH \\ \\0(nrH)\\ 

+ \\K{TH)\\\\k{fHr-0{ntH)\\ (36) 

< e(r/3) + (1 + 5(rP))(n£i(t f3) + z>(f/3)). 

Since, as it can be noticed in Table [lj 5(6) — e(9) <C 1, then we can take simply 

\\K(t H)K(t H) n - Q(tH) || < e(r j3) + ny(f f3) + z>(f/3). (37) 


It is worth remarking that such an approximation will require 2(nrh+m) + l real 
matrix-vector products of the form Hv, and thus is equivalent in complexity to 
the application of a (Chebyshev or Taylor) polynomial approximation of degree 
nrh + m. 


3.3 Flow of the algorithm 

Once a set of symplectic splitting methods constructed for providing approxi¬ 
mations under different conditions are available (methods collected in Table [T]) 
we still have to design a strategy to select the most appropriate scheme and 
step-size to carry out the numerical integration in time with the desired accu¬ 
racy and a as small as possible computational cost. 

The user has to provide the values for T m ; n and F max , a subprogram to 
compute the product Hv for a given real vector v, the final integration time t 
and the desired error tolerance tol. The procedure then implements the shifting 
<©, computes the value of (5 and determines the normalized Hamiltonian H. 

Next, the algorithm determines the most efficient method (or composition of 
methods) among the list of available schemes which provides the desired result: 
it chooses the cheapest method with error bounds below such tolerance and, if 
several methods with the same computational cost (same value of m) satisfy 
this condition, the algorithm chooses the scheme with the smallest error bound. 
This can be achieved if one starts the search from the methods with the smallest 
value of m and, for each value of m, proceeds by decreasing accuracy, i.e. by 
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±v±m 

m 

6» = 

/37m ax 

y*/m 

e(0) 

H(6) 

v(6) 

5(6) 

iw 10 

10 

5 

0.63 

3.6 x 10" 8 

8.7 x 10 n 

9.8 x 10" 8 

3.6 x 10 -8 

M (a9) 

iw 10 

10 

9 

0.94 

3.4 x 10" 5 

2.9 x 10" 5 

1.1 x 10“ 5 

6.0 x 10" 6 

M (°.e) 

iw 20 

20 

12 

0.79 

1.6 x 10 -13 

1.4 x 10 -13 

5.8 x 10 -14 

2.5 X 10“ 14 

ilJ 20 

20 

20 

1.1 

4.1 x 10" 7 

1.8 x 10" 8 

4.8 x 10“ 7 

4.0 x 10" 7 

M (a75) 

Jvi 30 

30 

22.5 

0.84 

8.1 x 10“ 15 

3.3 x 10“ 16 

1.5 x 10~ 14 

7.9 x 10- 15 

ivi 30 

30 

30 

1.0 

4.1 x 10 1 ° 

1.9 x 10 1 ° 

3.1 x 10~ 10 

2.6 x 10 1 ° 

iw 30 

30 

39 

1.36 

2.3 x 10" 5 

5.2 x 10“ 6 

2.2 x 10“ 5 

2.0 x 10" 5 

ilJ 40 

40 

40 

1.1 

1.8 x 10 -12 

4.9 x 10 -14 

2.4 x 10 -12 

1.8 X 10“ 12 

iw 40 

40 

48 

1.26 

2.1 x 10" 8 

2.1 x 10" 8 

5.3 x 10" 10 

4.7 x 10“ 10 

iw 40 

40 

56 

1.48 

1.48 x 10 -5 

4.0 x 10 -6 

1.7 x 10 -5 

1.7 x 10 -5 

iw 50 

50 

50 

1.07 

4.5 x 10 -15 

4.5 x 10 -15 

2.0 x 10 -17 

1.8 x 10~ 17 

1V1 50 

50 

55 

1.13 

4.5 x 10“ 13 

4.2 x 10“ 13 

4.1 x 10“ 14 

3.5 x 10” 14 

iw 50 

50 

60 

1.26 

5.4 x 10" 11 

2.7 x 10 -11 

3.8 x 10 -11 

3.4 x 10~ n 

, y.(1.3)a 
JU 50 

50 

65 

1.32 

1.2 x 10" 8 

1.2 x 10" 8 

8.3 x 10~ 10 

7.6 x 10 1 ° 

J\//( L3 b 

m 50 

50 

65 

1.32 

5.9 x 10“ 7 

9.5 x 10 -11 

6.1 x 10“ 7 

5.9 x 10“ 7 

iw 60 

60 

66 

1.15 

7.2 x 10" 15 

7.2 x 10“ 15 

2.6 x 10~ 17 

2.2 x 10” 17 

^(1.2)a 

ili 60 

60 

72 

1.3 

1.5 x 10 -12 

1.1 x 10 -12 

8.3 x 10 -13 

7.5 x 10“ 13 

^(1.2)fe 

ili 60 

60 

72 

1.26 

4.2 x 10" 11 

6.5 x 10" 14 

4.6 x 10" 11 

4.2 x 10“ n 

iw 60 

60 

78 

1.36 

1.2 x 10" 9 

7.8 x 10 -11 

1.2 x 10 -9 

1.2 x 10 -9 

j\ / j( 1 - 4 ) a 

1VI 60 

60 

84 

1.41 

8.4 x 10“ 8 

2.4 x 10“ 8 

7.4 x 10" 8 

7.1 x 10" 8 

^(1.4)b 

ili 60 

60 

84 

1.46 

2.9 X 10 -6 

3.7 x 10 -9 

2.9 x 10 -6 

2.9 x 10 -6 

Strang 

1 

1 

2 

1.8 x 10” 1 

4.7 x 10“ 2 

1.5 x 10“ 4 

1.3 x 10” 1 

Strang 

1 

1.4 

2 

5.1 x lO^ 1 

1.5 x 10 1 

4.0 x 10 1 

4.0 x 10 1 

Strang 

1 

1.9 

2 

1.34862 

0.606472 

2.4894 

1.1746 


Table 1: Relevant parameters of several symplectic splitting methods especially 
designed to integrate the semi-discretized Schrodinger equation using a time 
step t = t/n with a maximum value r max . Here y* stands for the stability 
threshold and e(6), y(0), v(0), and 5(6) (for 6 = /lr max ) are the coefficients 
(appearing in the error estimates obtained in Subsection 3.2) given in (29), 
), and (|35|) respectively. 
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increasing the value of 0 = /3r max . For a given value of t.(3 and tol the algorithm 
checks for each method if t/3 < /3r max and, if this condition is satisfied, then it 
examines if e(0) <tol. This procedure corresponds to the sequence of methods 
collected in Table [l] from top to bottom. 

If none of the methods from the table satisfy both conditions for t/3 and tol, 
then the time integration is split, i.e. t/3 is divided and a composition of one or 
several methods is used instead. Due to the high performance of the methods 
with the largest number of stages (in this case 60) the algorithm examines the 
cost of n steps for the six 60-stage methods where n = Floor [f/3/r max /3] and 
the last step is carried using one method from the list of methods. It chooses 
the cheapest methods with the smaller error bound among the composition of 
methods which provide the desired accuracy. 

In this way, if we denote by K$) the matrix associated to method Mm\ 
then the resulting splitting method corresponds to the composition 

(38) 


where the algorithm chooses the methods (labelled by 71,72, m), the time steps, 
r, f, and the value of n±, where n\ = 0 if the method uses just one step. If 
ri] > 0 the error bound is given by (37) while for m = 0 the error bound is just 
e(r/3). 

This strategy has been implemented as a Fortran code which is freely avail¬ 
able for download at the website [20], together with some notes and examples 
illustrating the whole procedure. 

In order to compare the efficiency of the resulting algorithm with the poly¬ 
nomial approximations based in Taylor and Chebyshev, with the error estimates 
collected in Table [I] we have represented in Figure [l] the computational work 
(equivalent to a polynomial approximation of degree m) required for different 
tolerances and values of /3t. Diamonds, squares and circles correspond to the 
error tolerances 10 -4 , 2 • 10“ 7 and ICC 11 , respectively, obtained with one or 
several steps of schemes in Table [l] Notice that our algorithm based on sym- 
plectic splitting methods provide better accuracy with a considerably reduced 
computational effort. 


4 Numerical examples 

Next we apply the algorithm based on symplectic splitting methods presented in 
section [3] to two different examples and compare its main features with Cheby¬ 
shev and Taylor polynomial approximations. For the first example, previously 
considered in m to illustrate Chebyshev and Lanczos approximations, we pro¬ 
vide in addition the codes we have produced to generate the results and figures 
collected here. These can be found at f20|. The second example illustrates the 
performance of the methods on a one-dimensional Schrodinger equation with a 
smooth potential. 
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Example 1. The problem consists in computing u(t) = exp(—itH)uo with 
uq E C n a unitary random vector and the tridiagonal matrix 


/ 2 - 1 
-1 2 -1 


H = - 
2 


V 


-1 2 -1 
-1 2 / 


E 


tNxN 


(39) 


The eigenvalues of H verify 0 < E^ < 2 for all k, so that we can take E m ; n = 0, 
i?max = 2, and thus a = /3 = 1 in Q. In consequence, the problem reduces to 
approximate 

where H = H — I . (40) 


e- iat e- ifHH UOj 


We take N = 10000 for the numerical experiments, but the results are largely 
independent of N (this is so even for the simplest, scalar case IV = 1). 

Both Chebyshev and Taylor methods have been implemented in such a way 
that only real valued matrix-vector products are used (we always separate into 
the real and imaginary parts, i.e. Hu = H(q + ip) = Hq + iHp )), so that 
Chebyshev requires to store only 7 real vectors instead of 4 complex vectors. 

We take as final time t = 20 and measure the error in energy, the error in 
the preservation of unitarity and the tolerance for different values of m, the 
degree of the corresponding polynomials. The results are shown in Figure [2] 
with the following notation: dashed lines for the relative error in energy, solid 
lines for the error in unitarity, and dotted lines for the theoretical error bounds 
of the approximate solutions. 

From the figure it is clear that the theoretical error bounds for the Taylor 
method are quite accurate for this example (since the bounds for E m ; n and 
E r nax are sharp) and that for the effective time-step t/3 considered, the error 
is exceedingly large for m below reaching the super linear convergence regime. 
This is not the case for the Chebyshev method (notice that the estimate © 
is valid only for m > t/3) since the coefficients c^- of the polynomial (18) do not 
grow as much as in Taylor. We also depict the results achieved by the first two 
splitting methods with r max /3 > t/3 = 20, and For these schemes 

the corresponding relative error in energy is represented by filled squares, the 
error in unitarity by filled circles and the error bounds by crosses. 

The relative performance of different numerical integrators is usually tested 
by measuring the error of the methods versus their computational cost. How¬ 
ever, the splitting methods we are considering in this work are designed to 
achieve a given tolerance, whereas their computational cost is determined through 
the error bound estimate. For this reason, we believe it is more appropriate 
to measure the cost of the methods for different values of the tolerance. In 
particular, we take tol = 10~ fc , k = 1,2,..., 12 and final integration times 
t = 20, 50,100, 200, 500,1000. Figure [3] shows the results obtained with Cheby¬ 
shev (line with squares) and the algorithm based on splitting schemes (line with 
circles) as a function of m. Even when high accuracy is required over long inte¬ 
gration times (the most advantageous situation for Chebyshev approximations), 
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0 10 20 30 40 50 60 70 


m 


Figure 2: Different approximations to e~ lt ^ H uo, with H given in (39)-(40), uq 
a random vector, f3 = 1 and t/3 = 20 versus the degree of the polynomials, 
m. The figure shows the relative error in energy (dashed lines), the error in 
unitarity (solid lines) and error bounds (dotted lines) for Chebyshev and Taylor 
methods. The results for the first two splitting methods with /3r max > /3t = 20, 
M 2 ( J } and are also shown: relative error in energy (filled squares), error 

in unitarity (filled squares) and error bounds (crosses). 


the new algorithm requires a smaller value of m and therefore less computa¬ 
tional effort. Notice how the algorithm selects the value of m to achieve the 
desired tolerance. 

Figure [4] shows the corresponding results for the relative error in energy 
versus m for the same example. Similar results are obtained for the error 
in unitarity or the two-norm error for which the error bounds apply (in this 
case one should compute numerically the exact solution and compare with the 
approximations obtained for each value of tol). 

Example 2 (Poschl—Teller potential). To illustrate how the methods work 
on a more realistic case, we consider the well known one-dimensional Poschl- 
Teller potential, which is an anharmonic quantum potential 

V( X ) = -?L- A(A ~ 1} 

2 fi cosh 2 (ax) ’ 

with a > 0, A > 1. It has been frequently used in polyatomic molecular simu¬ 
lation and is also of interest in supersymmetry, group symmetry, the study of 
solitons, etc. in is m- The parameter A gives the depth of the well, whereas a 
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t=20 t=50 t=100 





t=200 


t=500 


t=1000 





m 


m 


m 


Figure 3: Degree m of the polynomials to achieve tolerances tol= 10~ fc , k = 
1,2,..., 12 for different values of /3t (/3 = 1 for this problem) as determined 
by the error bound formulas using the Chebyshev method (squares) and the 
algorithm based on splitting methods (circles). 

is related to the range of the potential. The energies are 
a 2 

Ek = -(A — 1 — k) 2 , with 0 < k < A — 1. 

2 fj, 

We take the following values for the parameters (in atomic units, a.u.): 
reduced mass /i = 1745 a.u., a = 2, A = 24.5 (leading to 24 bounded states), 
and x E [—5,5]. Moreover, to apply a pseudo spectral space discretization we 
assume periodicity of the potential in this range. The resulting V(x) is thus 
continuous and very close to differentiable for all x E M. Table [2] collects the 
bounds to the spectral radius (obtained according to 0 ) and the corresponding 
shifting for the Poschl-Teller potential when the space interval x E [—5, 5] is 
split into N parts and for different values of N. Notice how sensibly E max 
depends on the space discretization. 

We take as initial condition a Gaussian function, ip(x, 0) = (re - ^ , where 
a is a normalizing constant, so the function and all its derivatives of practical 
interest vanish up to round off accuracy at the boundaries. The initial condi¬ 
tions contain part of the continuous spectrum, but this fact is largely irrelevant 
due to the smoothness of the periodic potential and wave function. 

Suppose that one is interested in solving the corresponding semi discretized 
problem in time with the following requirements: 

(I) N = 128, t = 157r, tol = 10~ 9 . In this case t/3 = 26.4648. 

(II) N = 512, t = 407r, tol = 10~ 6 . Now t/3 = 507.254. 
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Figure 4: Same as Figure [3] but replacing the value of the tolerance tol by the 
relative error in energy. 


We have to determine first, of course, the degree rn of the polynomial from 
the corresponding error bounds (for Taylor the time interval is divided by two 
in (I) and by 36 in (II) to avoid exceedingly large round off errors). Table [3] 
shows the number of matrix-vector products used by each method (in bold) and 
the 2-norm error for each method (compared with the exact solution obtained 
numerically with very high accuracy). In the first case our algorithm makes 
the computations in a single step using M while in the second case it uses 6 
steps of the scheme Adg' 4 ' 1 " followed by one step of m[q' 5 \ i.e. the composition 


(38) is now 


K 


(0.5), 
10 


(r/3) [Kg A)a (t 


with f = 84//3 and r = 407r — 6 t, and for a total of 370 products. Again, the 
algorithm based on symplectic splitting methods is able to produce results with 
the required accuracy with less computational effort. 
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N 

-S'min 

max 

a 

P 

64 

-0.65988 

0.11583 

-0.27202 

0.38785 

128 

-0.65988 

0.46333 

-0.098275 

0.5616 

256 

-0.65988 

1.8533 

0.59672 

1.2566 

512 

-0.65988 

7.4133 

3.3767 

4.0366 

1024 

-0.65988 

29.653 

14.496 

15.156 


Table 2: Bounds to the spectral radius and shifting for the Poschl-Teller po¬ 
tential with the parameters considered in the text, when the space interval 
x £ [—5, 5] is split into N parts. 



Taylor 

Chebyshev 

Symplectic 

tp = 26.4648 
tol = 10“ 9 

104 

3.4 x 10“ 12 

51 

3.7 x 10~ 12 

30 

4.2 x 10~ n 

tp = 507.254 
tol = 10 -6 

1836 

2.5 x 10~ 8 

587 

3.4 x 10” 15 

370 

4.4 x 10 -9 


Table 3: Number of matrix-vector products (in bold) and actual errors given by 
the Taylor, Chebyshev and symplectic methods for different t (3 and tolerances 
tol. 
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Appendix: Construction of methods 


We next describe the algorithm used to determine the coefficients (26) of length 
2m + 1 for given m and 9 G (0, 2m). 

Since all the error estimates in Subsection [372] depend exclusively on the even 
polynomial (of degree 2m) C(y) and the odd polynomial (of degree 2m+l) S(y) 
given in ( |30[ ), we first try to determine an appropriate pair of such polynomials 
satisfying the necessary conditions C(0) = 1 and C(x) 2 + S(x ) 2 — 1 > 0 (for 
all x G M). Such pair of polynomials is uniquely determined by a polynomial 
P(y) = C(y) + S(y) of degree 2m + 1 satisfying 


P( 0) = 1, l -(P(y? + P{-y) 2 )- 1>0. 


(41) 


Once an appropriate polynomial P(y) = C(y) + S(y) satisfying (41) is cho¬ 
sen, there is only a finite number of corresponding sequences (26), which can 


be effectively determined [2]. Since all of them share the same error estimates, 
we choose among them a sequence that minimizes 


ra +1 7 

£ M + £ \ h 


’j\- 


3 = 1 


3 = 1 


We next focus on the effective construction of the polynomial P(y) = C{y) + 
S(y) of degree 2m + 1. 

On the one hand, in order that the expression \JC\y) 2 + S(y) 2 — 1 featuring 
in the error estimate (29) be small in the interval y G [— 9 , 9\, 


sup 

-0<y<6 


cos (y + e(y)) + sin(y + e(y)) - P(y)\ 


(42) 


should be small for some real valued function e(y). On the other hand, mini¬ 
mizing 

V (C(y) - cos (y)) 2 + (S(y) - sin (y)) 2 ) 

in the interval y G [— 9,9\ is, provided that ( f42| is small enough, essentially 
equivalent to minimizing 

sup |e(y)|. (43) 

~e< y <o 

To reduce the complexity of the final algorithm for determining the polynomial 
P(y), we will try to minimize instead an alternative norm of e{y) that we 
introduce next. First observe that if 


e(y) = e 0 + ^2 % T i ( yI 6 ) ( 44 ) 

3> 1 


is the Chebyshev series expansion of the function e(y), then 


sup 

-e<y<6 


e{y)\ < 

3> 0 


-3 I' 


(45) 
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This suggests that the right hand side of (45) may be a good alternative to the 


supremum norm for sufficiently smooth functions e(y). For practical consider¬ 
ations, we will minimize instead the following alternative norm of the function 

e(y) 

(46) 



e g = 


Now, to determine the polynomial P(y) = C(y) + S(y) of degree 2m + 1, 
we consider, for a given odd integer l such that m + 1 < l < 2m, a given set 
of nodes yi,--.,yi symmetrically placed in the interval [—9, 9\, and a given odd 
polynomial e(y) of degree 1—2, the polynomial P(y) of degree 21 — 1 interpolating 
in the Hermite sense the function cos (y + e(y)) + sin (y + e(y)) for the nodes 
yi,... ,yi. In particular, this implies that P(0) = 1 and 

C(yf + S(y) 2 - 1 = \(P(yf + P(-y) 2 ) - 1 = V(y)W(y) 2 (47) 

where W(y) = (y — y\) • • ■ (y — yi), and V(y) is an even polynomial of degree 


4m — 21 + 2. Thus, P(y) satisfies the necessary condition (41) if and only if 
V(y) > 0 for all y. 


Notice that the interpolation error (42) admits an upper bound of the form 


sup | cos (y + e{y)) + sin(y + e(y)) - P(y)\ < 

-6<y<0 


-J— sup W(y) 2 , (48) 

[zt)\ -e< y <e 


where rj > 0 is an upper bound of the (absolute value of) the (21)th derivative 
of the function cos (y + e(y)) + sin (y + e(y)) in the interval y £ [—6,0\. 

For a prescribed set of nodes y±,... ,yi, we restrict the choice of the odd 
polynomial e(y) (of degree 1 — 2 ) so that the Hermite interpolating polynomial 
P(y) is of degree 2m + 1 (which introduces 2(1 — m) — 2 non-linear constraints 
on the non-zero coefficients e±, <33 ,..., e/ of the polynomial e(y) given by ([44])), 
and determine e(y) by minimizing the norm ||e||0 for that restricted set of odd 
polynomials e(y) of degree l — 2. This produces a polynomial P(y) for each 
choice of the set of nodes y\,... ,y[. It then remains to choose, for a prescribed 
positive odd integer l, such a set of nodes y±,... ,yi. 


The error estimate (48) suggests that a good choice for the interpolating 
nodes {yi, ■ ■ ■ ,yi} may be given by the zeros of the Chebyshev polynomial 
Ti(y/6 ) of degree l, which corresponds to minimizing the supremum norm (in the 
interval [—9, 6 ]) of the polynomial W(y). Notice that minimizing the alternative 
norm ||W||0 also gives rise to the same set of nodes. It then only remains, for 
given odd positive number 2m +1 and for given 9 > 0, to determine the number 
l of interpolating nodes, that should satisfy m + 1 < l < 2m. If l is too close to 
2m, then, very few degrees of freedom are left to minimize ||e||0, and if l is too 


close to m + 1, then the Hermite interpolating error (42) is too large, causing 
the norm of the function C(y) 2 + S(y ) 2 — 1 not being small enough, in addition 


to V(y) in (47) typically not being positive. We thus proceed by determining 
P(y) = C(y) + S(y) for different values of l close to (3m + 3)/2, and choosing, 
among those satisfying V(y) > 0, one having the best error coefficient e(9) 


defined in (29) 
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Unfortunately, choosing the interpolating nodes {yi,..., yi} as the zeros of 
the Chebyshev polynomial T)(y/0) of degree l typically results in a polynomial 
P(y) = C(y) + S(y) that does not satisfy the stability condition 


\C(y)\ < 1, ye [-9,0], 


(49) 


so that the error coefficients y,{Q),v{6) are not well defined, and thus the re¬ 
sulting splitting method cannot be reliably used in a step-by-step manner for 
large values of t/3. In order to produce splitting methods satisfying that sta¬ 
bility condition for given 9, we proceed iteratively to choose the interpolat¬ 
ing nodes {yi, ■ ■ ■ ,yi} and the corresponding polynomial P(y) as follows: As 
a first approximation, we require the set of nodes {y±,... ,yi} to contain the 
set {j tt : j e Z, |j7r| < 9} and determine the remaining nodes by mini¬ 
mizing the norm ||W||e of W(y) = (y — y\) ■ ■ ■ (y — yi). Once the polynomial 
P(y) = C(y) + S(y) is determined for that set of nodes {yi,..., yi}, we com¬ 
pute the set of zeros of C'(y )) = 0 that are included in the interval [—9,9] (that 
are typically close to {jir : j e Z, |j7r| < 0}), and determine the remaining 
nodes by minimizing the norm || W||# of W(y) = (y — yi) • • • (y — yi). Successive 
iteration of this process gives a sequence of polynomials P(y) = C(y) + S(y) 
that converge to a polynomial satisfying the stability condition (49j). 

As an example, we have obtained the method Mg q’ 4 -“ in Table TTby following 
this procedure for m = 60, 9 = 84, and l = 97, which has produced a splitting 
methods with sequence of coefficients (26) plotted in Figure [5] 


0.10 



Figure 5: Graphical representation of sequence (ai, b\, ct2, 62, • • •, ^60, &6O1 &6i) 
of method Mgg' 4 ^' in Table [l] obtained with 9 = 84 and l = 97. 
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